Quantum Monte Carlo simulations of confined bosonic atoms in optical lattices 
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We study properties of ultra-cold bosonic atoms in one, two and three dimensional optical lat- 
tices by large scale quantum Monte Carlo simulations of the Bose Hubbard model in parabolic 
confinement potentials. Local phase structures of the atoms are shown to be accessible via a well 
defined local compressibility, quantifying a global response of the system to a local perturbation. 
An indicator for the presence of extended Mott plateaux is shown to stem from the shape of the 
coherent component of the momentum distribution function, amenable to experimental detection. 
Additional fine structures in the momentum distribution are found to appear unrelated to the local 
phase structure, disproving previous claims. We discuss limitations of local potential approxima- 
tions for confined Bose gases, and the absence of quantum criticality and critical slowing down in 
parabolic confinement potentials, thus accounting for the fast dynamics in establishing phase coher- 
ence in current experiments. In contrast, we find that flat conflnement potentials allow quantum 
critical behavior to be observed already on moderately sized optical lattices. Our results furthermore 
demonstrate, that the experimental detection of the Mott transition would be significantly eased in 
flat confinement potentials. 

PACS numbers: 03.75.Hh,03.75.Lm,05.30. Jp 



I. INTRODUCTION 

Experiments on ultra cold atomic gases in optical lat- 
tices 0, provide the unique opportunity to directly 
compare theoretical studies of strongly correlated many 
body quantum lattice models with near-perfect experi- 
mental realizations of those models . Unlike for other 
strongly correlated systems, where often the models sim- 
ulated numerically are prototype toy models for describ- 
ing the properties of such materials, the same models are 
here realistic models of confined cold atoms, allowing for 
quantitative comparisons. 

So far, experiments have focused on bosonic atoms 1, 
0, which are easier to cool than fermions. A quantitative 
understanding of these bosonic systems will, in the fu- 
ture, allow to correctly interpret measurements on atomic 
fermion gases, which could be used as analog quantum 
simulators for strongly correlated fermionic systems Q . 

Here we present results of an extensive set of quantum 
Monte Carlo simulations performed on one dimensional 
(ID), two dimensional (2D) and three dimensional (3D) 
systems of confined bosonic atoms in optical lattices, ex- 
tending previous work in ID |^ and 3D [gl. Parts of 
our results were already presented elsewhere Q , here we 
provide a more detailed discussion and new results on 
confined Bose gases. Details on the many body quan- 
tum lattice model used for simulating confined bosons, 
as well as the employed quantum Monte Carlo method 
are discussed in the following section. 

In the first part of the paper, Sec^I we focus on simu- 
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FIG. 1: Spatial dependence of the local density for two di- 
mensional confined atomic boson systems: a) in the superfluid 
phase; b) at stronger lattice depths a central Mott insulat- 
ing region is formed. The simulations were performed for a 
parabolic trap with curvature V/U = 0.002, at ^/U = 0.37, 
for U/t = 6.7 (a), and U /t = 25 (b). 



lations of bosons in 2D confinement potentials, and study 
various aspects of confined bosons in a realistic setup. In 
particular, we (i) identify a local probe for the state at a 
given location inside the inhomogeneous trap, (ii) present 
a detailed analysis concerning the limitations of local po- 
tential approximations, and (iii) discuss the nature of 
spatial correlations inside the trap. 

Next, in Sec. IIVI we introduce an effective model de- 
scribing the inhomogeneous Bose gas inside the confine- 
ment potential. This allows us to study the finite size 
effects arising from the limited extent of the Bose gas in 
the trap. Using this effective model, we discuss the ab- 
sence of quantum critical behavior for bosons trapped in 
confinement potentials. 

In a third part. Sec. El we aim at identifying experi- 
mentally accessible quantities that signal restructurings 
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in the density distribution of the confined Bose gases. For 
this purpose, we monitor the evolution of the momentum 
distribution function upon increasing the lattice depth 
of the optical lattice, similar to the experimental proce- 
dure l] . We find that a method put forward recently (Qj 
for identifying the phase structure of the confined Bose 
gas is not appropriate. However, by monitoring the evo- 
lution of both the coherence fraction and the full width 
at half maximum of the coherent part of the momentum 
distribution function we obtain clear signals for restruc- 
turings of the density distributions. Finally, we summa- 
rize our observations in Sec. IVII along with implications 
for the experiments in confined atom systems. 



II. MODEL AND NUMERICAL TECHNIQUES 

Cold confined atomic Bose gases in an optical lattice 
are direct experimental realizations of the inhomogeneous 
Bose Hubbard Hamiltonian 3] 



(1) 



where 6j, (bi) denotes the creation (destruction) operator 
for bosons at lattice site i, located a distance ri from the 
center of the trap, and rii — b\bi the local density opera- 
tor. The nearest neighbor hopping integral t, the onsite 
Hubbard repulsion U , and the curvature of the parabolic 
confinement potential ^ > are tunable parameters in 
experiments. 

After evaporative cooling of the atoms, the experi- 
ments are, to a good approximation, performed at con- 
stant particle number. Similarly, the chemical potential 
/i allows us to control the filling of the trap in the numer- 
ical simulations. For the discussions below we introduce 
a local chemical potential 
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which decreases upon moving away from the trap center. 

Quantum Monte Carlo (QMC) simulations of the 
Hamiltonian Eq. |^ were performed using the Stochastic 
Series Ex pan sion method ^8, 9J, with directed loop up- 
dates HJI ■ This algorithm requires a cutoff N^^y, 
on the local site occupation. Performing simulations at 
mean densities (rzi) < 1, a cutoff iVmax = 2 or 3 can be 
chosen without introducing significant errors. The tem- 
perature used in the QMC calculations was chosen low 
enough to be essentially in the ground state. 

All lengths are set in units of the lattice constant a, 
and the simulation box is taken large enough, to ensure 
that outside its boundary the local density rii vanishes 
due to depopulation by large negative values of . For 
the values chosen in our simulations we needed to keep 
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up to 500-site chains in ID, 50 x 50 square lattices in 2D 
and up to 16^-site cubes in 3D. 

To distinguish the local phases in the inhomogeneous 
system and to probe for quantum criticality we define the 
local compressibility at site i, 

(3) 

by the response of the system's density, n, to a local 
chemical potential change at site i. Here, (3 = l/T de- 
notes the inverse temperature. Another way of probing 
for local properties is to measure local density fluctua- 
tions These are related, by 



1 / dui 
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to a local density response. In the homogenous case, the 
global compressibility 
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is equal to the local comressibility nf"^^^ defined in Eq. Q 
at each lattice site i, but is in general not given by /3Ai. 
This is due to correlations between different lattice sites 
in the superfluid regime. Only in the Mott-insulating 
phases of the homogenous system, where these correla- 
tions are absent, do all quantities become zero. We will 
show in Sec. lHI Bl that the local compressibility, as 
defined in Eq. Q is indeed able to characterize the local 
state near a given site even in the inhomogenous case. 

In addition to measuring local quantities, we measure 
the momentum distribution function. 
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where N denotes the total number of particles within the 
system. This way the momentum distribution is normal- 
ized to unity, and the coherence fraction given by the 
height of the coherence peak, n(k = 0). 

The momentum distribution function n(k) is directly 
accessible in experimental studies of confined atomic 
gases. For negligible interparticle interactions during bal- 
listic expansion of the atomic cloud, it essentially maps 
onto the interference pattern of the resulting absorption 
imagines [T^ . Due to the tight binding approximation 
of the Bose-Hubbard model, the numerical momentum 
distribution functions are periodic in the extended zone 
scheme. A finite extent of the Wannier functions on each 
site in an optical lattice leads to a form factor in the mo- 
mentum distribution resulting in unequal heights of 
the higher order Bragg peaks, as observed in the experi- 
mental interference patterns yj. 

To further quantify spatial correlations within the trap, 
we measure the one particle density matrix, i.e. the cor- 
relation function 



gi^,J) = {blb,), 

between bosons on lattice sites i and j. 
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III. SIMULATIONS OF REALISTIC 2D TRAPS 



A. Phase coexistence in trapped systems 

Trapped inhomogeneous systems can show coexistence 
of both supcrfluid and Mott insulating regions, as has 
been clearly identified in experiments j in mean field 
investigations numercial renormalization group 
and Gutzwiller ansatz calculations [l5j and in numerical 
simulations in ID II III "l^ and 3D 6J. Simulat ing large 
realistic systems in 2D, we can clearly observe this coexis- 
tence. For large hoppings t the total system is superfluid 

(Fig.nji). 

As the hopping amphtude t is decreased, a Mott insu- 
lating plateau with integer density (here (rii) — 1) forms 
in the center of the trap, as seen in Fig. ^p. A super- 
fluid ring-like region with a non-uniform particle density 
surrounds the central plateau. The typical width of this 
superfluid ring is 6 — 8 times the lattice unit, for the 
parameters used in our simulations. Far away from the 
center, the density is (rii) — which can also be viewed 
as a Mott plateau. 

The emergence of the Mott plateau and the shrink- 
ing of the superfluid region have been interpreted as a 
quantum phase transition j|. We prefer to view it as a 
crossover where the volume fraction of the Mott insulat- 
ing phase grows and that of the superfluid phase shrinks. 
After discussing the quantitative properties of the two 
phases and of the boundary region we will, in Sec. IIV CI 
show that no quantum critical behavior is observed in 
this system, which supports our interpretation of this 
phenomenon as a crossover instead of a phase transition. 



B. Local compressibility 

Fig. n is useful for a qualitative view, but we need 
a quantitative probe to distinguish the different phases. 
This probe is given by the compressibility. QMC results 
for the local compressibility within the trap of Fig. 
are shown in Fig. [3 The extent of the superfluid shell is 
clearly reflected by this local probe, which also shows that 
the compressibility is largest close to the outer regions of 
the superfluid shell. 

For a quantitative analysis we plot in Fig. |31the radial 
dependence of both the local density n{r) and the local 
compressibility K'°'^'*'(r) as a function of the distance r 
from the center of the trap. We observe two well defined 
Mott regions with local density n{r) = 1 and n{r) = 0. In 
these Mott plateaux, the local compressibility vanishes, 
whereas it is finite in the intervening superfluid ring with 
non-integral density < 7i(r) < 1. 

More precisely, we observe two well defined peaks in 
K^°'^^^{r), signaling an increase in the particle number 
fluctuations at the boundaries to the Mott regions. While 
these two peaks would be of the same height in hard-core 
boson models (due to particle- hole symmetry), they are 
asymmetric in the soft-core case. 
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FIG. 2: Spatial dependence of the local compressibility, k "'^^ , 
of bosons in a two dimensional parabolic trap with curvature 
V/U = 0.002, for fi/U = 0.37 and U/t = 25. A superfluid 
ring surrounding the n — 1 central Mott plateau is clearly 
resolved. 




FIG. 3: Radial dependence of the density, n, and local com- 
pressibility, of bosons in a two dimensional parabolic 
trap with curvature V/U = 0.002, for fj./U = 0.37 and 
U/t = 25, with a superfluid ring surrounding the n = 1 cen- 
tral Mott plateau. The inset shows the radial dependence of 
the local density fluctuations. A, which remain finite within 
the Mott plateau region. 



The inset of Fig. O shows the radial dependence of the 
local density fluctuations of Eq. A(?^ flrst used in 
simulations of ID conflned Bose systems [3. As opposed 
to K}°'^^^{r), this quantity peaks in the middle of the su- 
perfluid ring, and does not vanish in the central Mott 
plateau region. While density fluctuations are therefore 
most pronounced in the middle of in the superfluid, they 
remain flnite inside the Mott plateau, due to the sur- 
rounding superfluid. These results demonstrate that the 
local compressibility K^°'^'^^{r) is a better probe for the ex- 
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FIG. 4: Spatial dependence of the correlation function 
g{hj) = {b\bj) between bosons at site Vj — (12.3, 0) inside the 
superfluid shell, and all other sites i within a two dimensional 
parabolic trap with curvature V/U — 0.002, at /i/U — 0.37 
and U/t = 25. 



istence of superfluid or Mott regions in the system than 
the onsite response expressed by A(r). Moreover, the 
response of the total system to a local excitation should 
be easier to study experimentally than the local response 
to a local excitation (for example by changing the laser 
intensity at one specific point). 



C. Spatial correlations in the superfluid 

In order to gain a better understanding of the extended 
superfluid ring, we study the behavior of the one parti- 
cle density matrix, Eq. Q, in the trapped system. In 
Fig. 21 the spatial dependence of g{i,j) is shown for a 
fixed site well inside the superfiuid ring, and all other 
sites i in the parabolic trap. The rapid decay of the cor- 
relations towards both the central Mott plateau region 
and the outside clearly exhibits the ring-like structure of 
the coherent superfiuid. 

Inside the superfluid ring the correlation function de- 
cays less rapidly. In order to make a quantitative analysis 
of its long distance behavior, we plot in Fig.[51the depen- 
dence of the correlation function g{i,j) for site i,j along 
the ring on the distance d between the sites, as measured 
along the ring. See the left inset of Fig. |31 for an illus- 
tration. We flnd all data for distances d/a > 5 to follow 
very closely a flnite size exponential decay. 



g{d) = c + 6 cosh 



(8) 



with a correlation length ^/a = 21.6 ± 0.4, and a finite 
constant c = 0.033 ± 0.003. Thus, the superfiuid ring 
does not exhibit quasi- long ranged correlations, as might 
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FIG. 5: Dependence of the correlation function g(d) = {b\bj) 
on the distance d between sites i and j along the superfluid 
ring for sites within r/a = 12.3 ± 1.2 lattice units from the 
center of a two dimensional parabolic trap with curvature 
V/U = 0.002, at ^./U = 0.37 and U/t = 25. For d/a > 5 
the data fit well (x^ = 1.6) to an exponential decay, with 
a finite value of c = 0.033 ± 0.003, and correlation length 
^/a = 21.6 ± 0.4. The right inset shows the correlation at 
the largest distance, g{nr), as a function of the distance r 
from the trap center. The left inset illustrates the distance d 
between two sites i,j along a ring of radius r employed here. 



have been expected from the reduced dimensionality of 
the quasi-lD superfluid ring. Instead, we find the ring to 
be wide enough for long range order, i.e. 2D behavior, 
to persist inside the superfluid in spite of the presence of 
a central Mott plateau. 

A rapid decrease of the correlations, as one moves away 
from the middle of the superfluid, can be seen from the 
right inset of Fig. [S] which shows the radial dependence 
of the correlation function g{d) at the largest distance 
d = 7rr. The correlations are thus most pronounced near 
r « 12, which corresponds to the central region of the 
superfluid ring, and become suppressed as one moves to- 
wards either boundaries. Finally, inside the Mott plateau 
region, the correlations are merely short ranged. Even 
though the coherent part of the Bose gas is thus conflned 
to a ID ring-like structure, the spatial dependence of the 
correlations inside the ring clearly exhibit the underlying 
2D structure of the trapping potential. 

Our flndings are related to recent results in ID sys- 
tems 0, which show that the presence of a Mott 
insulating region does not change the long distance be- 
havior of the one particle density matrix in trapped boson 
systems. Similar analysis might also apply to 3D con- 
finement potentials, which would allow for extended 2D 
superfiuid shells surrounding a Mott insulating central 
region. 
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FIG. 6: Local density, n, and local compressibility, k}°'^^^, 
of bosons in a two dimensional parabolic trap with curvature 
V/U = 0.002, at fj./U = 0.37 and U/t = 25, as functions 
of fjf^^ . The results for from the density fluctuations, 

Eq.(0 (circles) agree with the numerical derivative dn/dff^ 
(triangles). 



FIG. 7: Density, n, and compressibility «: of the Bose- 
Hubbard model on the 2D square lattice for U/t = 25, as 
a function of /i. The data are taken from simulations on a 
24 X 24 sites square lattice. 



D. Local potential approximation 

The confinement potential V enters the Bose Hubbard 
Hamiltonian, Eq. by coupling to the boson density. It 
therefore corresponds to a local chemical potential shift, 
expressed by the local chemical potential fx'^^ , defined in 
Eq. ^ . In Fig. El we show the behavior of both the local 
density and compressibility as functions of this effective 
chemical potential, The observed data collapse in 

both quantities, as well as the smooth behavior of these 
curves for all points on the lattice, indicate that for the 
realistic parameters used for the simulation, a local po- 
tential approximation holds, i.e. the local density can be 
determined from the value of the local chemical potential. 
The validity of this approximation is further confirmed 
by observing that the local compressibility coincides per- 
fectly with the numerical derivative of the curve n{ij,°^), 
i.e. 

It is important to note, that n{n°^) is not a universal 
function, but depends on both the confinement geome- 
try as well as the trap curvature, V. Universal behavior 
is violated in particular near the boundary of the super- 
fluid, where n{ii'^^) does not behave as n{fi) for uniform 
systems, corresponding to = 0. For example, in the 
uniform 2D case a cusp develops when n approaches 1 , as 
seen in Fig.[7| due to quantum criticality. This cups is ab- 
sent in the confined system, where n approaches 1 rather 
smoothly (c.f. Fig. While a local potential approx- 
imation in terms of an effective local chemical potential 
holds, the value of n(/i°*^) thus has to be taken from simu- 
lations in a trap, and not from the uniform system on the 



underlying lattice. Approximative schemes such as those 
of Ref. 19], which assume that locally the system can 
be mapped onto an unconfined system at the local value 
of n'^^ , have to be applied with care. Being reasonable 
in the bulk of both the superfluid and the Mott plateau 
region, such approaches become unreliable near the in- 
terfaces between superfluid and Mott insulating regions, 
where differences between the confined and unconfined 
case are most pronounced. There QMC simulations are 
needed in order to obtain n{iJ,°^) for a given trap geom- 
etry. 



IV. THERMODYNAMIC LIMITS AND 
QUANTUM CRITICALITY 

In the previous section we found that in the boundary 
layer between Mott plateaux and superfluid regions, the 
confined system does not show cusps in the density pro- 
file like the uniform system does as a function of fj, near 
the quantum critical points. This is indication for the 
absence of quantum criticality in the confined system. It 
could be argued that this is simply due to finite size ef- 
fects, but here we show that the situation is indeed more 
subtle. For this purpose, we introduce an effective ladder 
model for bosons inside inhomogeneous potentials. 



A. Effective ladder model for the superfluid ring 

A parabolic potential imposed on a regular hyper-cubic 
lattice makes for an irregular form of the superfluid ring 
surrounding the Mott insulator (see Fig. . We thus 
flrst investigate whether this structural quenching is the 
reason for the apparent absence of quantum critical fea- 
tures by looking at a "smoother" lattice. Although it 
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FIG. 8: DifFerent lattice topologies and expected superfluid 
regions (denoted in grey) for a quadratic confining potential 
for such topologies: (a) square lattice, (b) polar lattice, (c) 
ladder lattice. 



cannot be realized in experiments, we consider a polar 
lattice, such as depicted in Fig. |S1 which is topologically 
equivalent to the ladder system in Fig.|Sl The rotational 
invariance of the polar lattice (translation invariance of 
the ladders along the leg direction) simplifies the detec- 
tion of critical behavior by reducing finite size effects. 

While the polar and ladder lattices are not topolog- 
ically equivalent to the square lattice, universality of 
quantum critical behavior ensures that we study the same 
effects. Comparing the lattices will be a further check 
for the validity of a local potential approximation for the 
confined system. 

The Hamiltonian of the ladder system of size Lx x Ly 

is 

= E E (4.^^+1.^ + + ^-c-) (10) 

1=1 j=i 

jj Lx Ly Ly Lx 

j=l i=l j = l 1=1 

where the confinement potential is now given by a chem- 
ical potential fi{j) which is constant for each leg j of 
the ladder. The boundary conditions are chosen to be 
periodic along the legs of the ladder, and open in the 
transverse direction. 

With these choices, the local density (riij) is indepen- 
dent of i due to translation invariance along a leg of the 
ladder, and we denote by n(j) the density of particles on 
the j-th leg. Similarly, the local compressibility K^°'''^^(j) 
is the same for all sites on the j-th leg. The results shown 
below have been obtained for a geometry with = 64, 
and Ly = 10. 

To simplify simulations, we linearize the quadratic po- 
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FIG. 9: Local density, n, and local compressibility, fi:'°'^^', for 
a two dimensional parabolic trap on a square lattice and for a 
ladder, as functions of ^'^^ . The parameters for the trapping 
potential are V/U — 0.002, at n/U — 0.37, and the param- 
eters of the ladder model are chosen as to cover the whole 
superfluid region (see text). The density profile within the 
trap and the ladder coincide almost perfectly. Small differ- 
ences are exhibited by the two different compressibility curves, 
which nevertheless share the same overall shape. 



tential, fi + Vj^ and set 

Kj) = ^^o+jA^l, (11) 

where Afi is the difference in chemical potential be- 
tween two neighboring legs. In our simulations we fix 
A^/U = 0.053 and make sure that the first leg is always 
in the {n) = Mott insulating phase and the last leg 
in the (n) = 1 insulating phase. Sweeping the value of 
the global /io allows us to obtain results for all values 
of the local chemical potentials, and to drive one of the 
legs across the transition from the superfluid to the Mott 
insulating phase. We note, that a similar construction is 
obtained for 3D traps, by modeling the superfluid shell as 
a set of coupled planes with a chemical potential gradient 
applied perpendicular to the planes. 



B. Results of the ladder model and comparison to 
the realistic trap 

In Fig. iniwe show the combined results of all simula- 
tions with different values of /io- We plot both the local 
density {n,{j)) and local compressibility K'°'^'''(j) for all 
legs and all values of /io as functions of the local chemi- 
cal potential. Note that all data from simulations using 
different values of /io collapse onto a single curve. Due to 
the equivalence of all sites on the ladder with the same 
value of the chemical potential, this data collapse is ex- 
pected in the ladder model, where changing /io by an 
amount of A/i corresponds to a shift in the index of each 
chain in the ladder. 
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Fig. 1^ also shows the results obtained for the realistic 
trap, as described in the previous section. Both results 
coincide almost perfectly, which is a strong justification 
for the validity of the effective model and the applicability 
of a local potential approximation. 

The small differences can be attributed to the different 
shape of the trapping potential which is parabolic for the 
realistic trap, but linearized in the ladder model. While 
differences in the densities are small, the compressibility, 
which is the derivative of the density, is more sensitive, 
marking the lack of universality. 

In particular, we find that the ladder model does not 
show cusp-like features in the boundary regions between 
the superfluid and the Mott plateaux, similar to the real- 
istic 2D trap. We thus conclude, that structural quench- 
ing in the rotational symmetric trap on the square lattice 
is not the reason for the apparent absence of quantum 
criticality since similar behavior is also found on the po- 
lar lattice, i.e. in the ladder model. We expect similar 
features also to be obtained for 3D traps using the effec- 
tive multilayer model. 



C. Quantum criticality 



ID traps compared to the uniform case 




FIG. 10: Local density, n, for bosons confined in one dimen- 
sional parabolic confinement potentials of different curvatures 
V, as a function of //'^ . For comparison, the behavior of the 
density as a function of ^ in the uniform ID system is also 
shown. The insets focus on the regions close to n = and 
n = 1, where diflFerences between the trapped case and the 
uniform case are most pronounced. 



Having shown that the ladder model captures impor- 
tant features of the realistic trap even quantitatively, we 
now discuss possible quantum criticality in confined sys- 
tems. 

Performing the thermodynamic limit in the uniform 
case does not change the model parameters apart from 
increasing the number of lattice sites. On the other hand, 
standard thermodynamic limit definitions for the con- 
fined case E] imply a decrease in the trapping po- 
tential's curvature V, such that 

TV -> oo, with N^/Vjt = const. (12) 

Such a procedure locally drives the system towards the 
uniform regime, since for lattice sites with a given value of 
the gradient in the confinement potential eventually 
becomes irrelevant on the scale set by the correlation 
length in their neighborhood. Therefore, the state of 
the uniform system at the value of is established 
there 113 . 

An example of this approach to the uniform system 
upon decreasing the trapping potential V is given in 
Fig. ^1 It shows the density n as a function of 
for two different ID traps. In the ID case, we can study 
larger systems than in 2D, allowing us to decrease V by a 
factor of 25, which requires increasing the linear system 
size by a factor of 5. While in both traps in Fig. llUlthe 
value of the local density follows the curve n(/i) of the 
uniform ID system, deviations are visible in the insets, 
which focus on the regions close to the Mott plateaux. In 
the case of the shallow trap, these deviations are clearly 
reduced, with more pronounced singularities developing 
at the critical point. 



The thermodynamic limiting procedure in Eq. H12|l cor- 
responds to the limit A/i — * in the gradient of the ef- 
fective ladder model. In this limit, we recover the Bose 
Hubbard model on the 2D square lattice, showing 2D 
quantum critical behavior. A finite gradient A/z in the 
ladder model restricts the correlation length perpendic- 
ular to the chains, and does not allow for 2D quantum 
critical behavior to be observed. 

However, in the opposite limit, A/i — > cxd, we recover 
a ID chain on the lowest leg, which will show ID criti- 
cal behavior. In the inset of Fig. the compressibility 
of the Bose Hubbard model is shown on a chain of 64 
sites. Already on such a short chain sharp peaks near 
the transition are visible, which diverge with increasing 
chain length p3 |. 

One might expect this ID quantum criticality to per- 
sist in the ladder model with a finite gradient Afi, ob- 
served upon changing the chemical potential hq to drive 
one of the legs across the phase transition. This is how- 
ever not the case: While the broad peaks in Fig. are 
weak remnants of ID quantum criticality, they do not 
diverge as the size of the ladder model is increased. This 
is seen from Fig. ^2 which shows the density and local 
compressibility as functions of for ladders of different 
lengths. 

Therefore, the effective ladder model as well as the re- 
alistic trapped system do not show quantum criticality. 
The quantum critical behavior that might have been ex- 
pected to occur in the boundary layer between the Mott 
insulator and the superfiuid is destroyed by the finite 
gradient in the effective chemical potential, and the cou- 
pling of this layer to the rest of the system. The absence 
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FIG. 11: Local density, n, and local compressibility, k "'^'^ , for 
the Bose Hubbard model on ladders with different lengths as 
a function of fi'^^ . The other parameters of the ladder model 
are chosen as in Fig. |^ The inset shows the compressibility 
for the uniform ID case as a function of /i for a chain with 64 
sites, for the same value of U /t as used in the ladder model. 
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FIG. 12: Momentum distribution function of bosons in a 
three dimensional parabolic trap with curvature V jlJ = 0.008, 
along the line (0, 0, 0)-(7r/a, 0, 0) in momentum space, for 
/i/f/ = -0.2, and U jt = 24, values taken from Ref. 0. The 
lower inset exhibits the presence of satellite peaks in this sys- 
tem, without a Mott plateau, as seen from the radial density 
distribution shown in the top inset. 



of quantum criticality in the realistic 2D parabolic trap 
is not only due to its finite size, but is in fact imposed 
by the inhomogeneity of the confinement potential, as 
we showed using the effective ladder description. The 
Mott transition observed in experiments on ultra-cold 
atoms should not be viewed as a quantum phase tran- 
sition, but instead as a crossover with changing volume 
fractions of the Mott plateaux and superfluid regions. 

In contrast, for flat confinement potentials, corre- 
sponding to y = 0, quantum critical behavior can be 
observed already on moderate system sizes. For exam- 
ple, the data shown in Fig. [7| for a 24 x 24 sites 2D sys- 
tem, and in the inset of Fig. ^] for a 64 sites chain, both 
clearly reflect the quantum criticality of the uniform 2D 
and ID cases. 



V. IDENTIFYING MOTT PLATEAU 
FORMATION IN TRAPS 

While measuring the spatial density distribution and 
the local compressibility in QMC simulations allows for 
the identification of the different regions inside the trap, 
such a local probe is not (yet) available experimentally. 

In order to identify changes in the density distribution 
inside the trap upon varying control parameters such as 
V ^ [/, or i, different strategies have been proposed. In 
particular, in Ref. @ it was claimed that the presence 
of a Mott plateau in the trap center is signaled by fine 
structures in the momentum distribution function, and 
that these are absent in traps without Mott plateaux. 

Two of the configurations of Ref. are shown in 
Figs. El and E| Here, we simulated 3D systems in a 
trap and used the parameters given in Ref. [6j, and plot 
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FIG. 13: Momentum distribution function of bosons in 
a three dimensional parabolic trap with curvature V jU = 
0.0122, along the line (0, 0, 0)-(7r/a, 0, 0) in momentum space, 
for /x/f/ = —0.3125, and U jt = 80, values taken from Ref. 
The corresponding radial density distribution is shown in the 
inset. 



the resulting momentum distribution along the line from 
k = (0,0,0) to (7r/a,0,0), along the k^-tcxSs of momen- 
tum space. 

While in Fig. 1131 where an extended Mott plateau is 
present, additional fine structures in n(k) are visible, 
such structures seem not to exist at first sight for the data 
shown in Fig. 1121 However, similar fine structures are also 
present there, even on a similar scale, as clearly seen from 
the bottom inset of Fig. ^1 which focuses on the tail of 
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the momentum distribution function. These fine struc- 
tures become invisible on the scale used for the main part 
of Fig. lEl due to the presence of a dominant coherence 
peak. We found in the simulations discussed below, that 
a simple correspondence as proposed in Ref. Q does not 
hold: The presence of fine structures is due to the finite 
extent of the superfluid within the trap, and emerges also 
without a Mott plateau being present in the trap center. 
Denoting the radial extent of the superfluid by R, these 
peaks appear for near integer multiplies of 27r/i?, if 
not masked by the incoherent background, or rendered 
almost invisible on the scale of the coherence peak. 

For a more systematic analysis of the experimentally 
accessible momentum distribution function we monitor 
its evolution upon varying experimentally accessible con- 
trol parameters of the system. We have performed sets of 
simulations, which mimic the experimental procedure of 
increasing the lattice depth of the optical lattice I'l, by 
reducing the hopping parameter t, while keeping all other 
parameters constant. We performed such scans for traps 
of different dimensionalities, and begin our discussion of 
the results for the 2D case. 
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FIG. 14: Density n as a function of U/t for bosons along a 
constant n/U = 0.37 scan, for both periodic (PBC) and open 
(OBC) boundary conditions on a 24x24 sites square lattice, 
OBC corresponding to an infinitely sharp trapping potential. 
The Mott transition in the uniform case is indicated by the 
dashed line. The inset locates this constant fj, scan within the 
phase diagram of the uniform system. 



A. Homogeneous and closed box 2D systems 

In order to better interpret data obtained for inhomo- 
geneous traps, we first consider the homogeneous case. 
To this end, we set the trap curvature V — 0, using 
periodic boundary conditions (PBC), and also consider 
a finite square lattice with open boundary conditions 
(OBC), representing an infinitely sharp trapping poten- 
tial, i.e. a closed box system. We then compare both the 
closed box and the parabolic trap to the uniform case. 

The phase diagram of the Bose Hubbard model on the 
square lattice in the parameter regime of interest is shown 
in the inset of Fig. E| While the overall phase structure 
is well described by mean field theory [23, |2J| , the maxi- 
mum extent of the first Mott lobe oiU/t = 23.3 is thereby 
overestimated by about 32% of the value obtained by a 
strong couphng expansion {U/t = 17.54) [^ilHl- Thus, 
in order to compare later the onset of Mott plateau for- 
mation in confined systems with the homogeneous case, 
we first need to determine the precise phase boundary 
within the region of interest. 

In the following we consider a path through the phase 
diagram along a line of constant ^/U = 0.37, indicated 
by the directed line in the inset of Fig.^] Moving along 
this line through the quantum phase transition, the sys- 
tem undergoes the generic transition 23] which is mean 
field in nature. The transition does not cross over to the 
special case of the commensurate transition, belonging to 
the 3D XY universality class [23j . 

In order to identify the Mott transition point along this 
scan we measured the stiffness, ps, which quantifies the 
response of the system to a twist in the boundary condi- 
tions along the real space directions. This quantity can 
be calculated in QMC from the boson winding number 
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FIG. 15: Scaling plot of the stifli'ness, ps, for the Bose Hub- 
bard model on the square lattice for periodic boundary con- 
ditions, while tuning along a constant p/U = 0.37 scan. The 
inset shows the stiffness, pa, for a 28 x 28 sites system as a 
function of U/t. 



fiuctuations 1261. At the quantum critical point finite size 
scaling theory [23j predicts it to scale in two dimensions 
as L~^, where z = 2 is the dynamical critical exponent 
for the generic transition, and L the linear system size of 
a. N = sites system. 

The inset of Fig. 1151 shows ps as a function of the in- 
verse hopping, U/t, for a system of linear size L = 28, in- 
dicating its increase upon entering the superfluid phase. 
The main part of Fig. ^| is a scaling plot of PsL^ vs. 
U/t, from which the critical point U/t = 16.25 ±0.1 is 
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FIG. 16: Momentum distribution function of bosons on a 
24x24 sites square lattice with periodic boundary conditions 
(PBC) for the commensurate momenta along the line (0, 0)- 
(tt/o, 0) in momentum space, for different values of U /t, while 
tuning t along the constant ^/U = 0.37 scan of Fig. 1141 The 
loss of coherence due to the Mott transition at U/t — 16.7 is 
reflected by the reduced coherence peak height, n(0, 0). 
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FIG. 17: Evolution of the coherence fraction, i.e. the height 
of the coherence peak, n(0,0), as a function of U/t while 
tuning t along the constant fi/U = 0.37 scan of Fig. 1141 for 
bosons on a 24 X 24 sites square lattice for both periodic 
(PBC) and open (OBC) boundary conditions. For OBC the 
full width at half maximum {FW H M) of the coherence peak 
is also shown, clearly signaling the onset of the Mott phase. 
The Mott transition point is indicated by the vertical dashed 
line. 



obtained with sufRcient precision for the purpose of this 
study. Note that corrections to scaling are visible already 
at L = 12, in agreement with a recent study, indicating 
that large systems sizes are needed for high precision de- 
terminations of the critical point Mean field the- 
ory [2^ predicts a value of U/t = 23.5 which is too large 
by more than 40%. Thus in 2D, corrections to mean field 
have to be accounted for already the uniform case. 

Having established the position of the Mott transition 
in the uniform case, we show in Fig. 1141 the evolution of 
the density, n, in homogeneous systems with PBC and 
OBC along the constant n scan shown in the inset of 
Fig. El Here, we show data on systems with 24 x 24 
lattice sites. In both cases the onset of the Mott insulat- 
ing regime is located close to the critical point obtained 
from finite size scaling. Thus already for this moderate 
system size, the Mott transition is located rather close to 
the value in the thermodynamic limit. 

Phase coherence is signaled by a pronounced coherence 
peak in the momentum distribution function, n(k — 0). 
This can be seen from Fig. 1161 which shows n(k) for the 
commensurate momenta along the /c^j-axis of momentum 
space for different values of iJ/f for a 24 x 24 sites system 
with PBC. The evolution of n(0) while tuning through 
the transition is shown in Fig. 1171 It clearly marks the 
loss of coherence in the Mott insulator. 

While for PBC a discrete set of commensurate mo- 
menta exists, for OBC states of arbitrary momenta can 
be occupied. In the QMC simulations, we measured 
the momentum distribution function on a mesh cover- 
ing 10 x L momenta in the first Brillouin zone along the 
/caj-axis. The resulting momentum distribution functions 
for different values of U/t are shown in Fig.^l Similar to 
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FIG. 18: Momentum distribution function of bosons on a 
24x24-site square lattice with open boundary conditions, cor- 
responding to an infinitely sharp trapping potential, along the 
line (0, 0)-(7r/a, 0) in momentum space, for different values of 
U/t, while tuning t along the constant fi/U — 0.37 scan of 
Fig. 1141 Within the superfluid regime satellite peaks appear, 
which diminish upon emergence of the Mott phase. 



the case of PBC, the loss of coherence upon entering the 
Mott phase is signaled by a reduction of the coherence 
fraction, n(0), as seen from Fig. El Also note the pro- 
nounced fine structures in n(k) for U/t = 10.0, deep in 
the superfluid regime, and the complete absence of such 
structures for larger values of U/t. 

Using a Ornstein-Zernike form for the coherent part 
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for the momentum distribution function, 



n(k) 



^(0) 

1 + 



(13) 



where ^ denotes the coherence length, we can obtain ^ 
from the fuU width at half maximum (FWHM) of the 
coherence peak, 



FWHM = 



(14) 



In the thermodynamic limit, the coherence length di- 
verges in the superfluid regime. On a finite system it 
is bounded from above by the linear system size L, and 
decreases to zero deep inside the Mott insulating regime. 
We thus expect an increase of the FWHM , from its min- 
imum value of 2 /L in the superfluid regime, upon driving 
the system through the Mott transition. This behavior 
can indeed be seen in Fig. 1171 which shows the FWHM 
as a function of U /t for the 24 x 24 sites closed box. The 
FWHM is at its lowest value of about 2/24 left of the 
Mott transition, and increases due to loss of coherence 
beyond this point. 



B. Parabolic traps in 2D 

After our analysis of the homogeneous system, we are 
now in position to discuss the evolution of a confined Bose 
gas in a 2D optical lattice while increasing the lattice 
depth, i.e. decreasing the hopping, t. In the following, 
we will take the system along a path of constant chemical 
potential, n/U = 0.37, in a parabolic trap of curvature 
V/U = 0.002. An illustration of the path taken in our 
simulations is shown in Fig. ^] 

As discussed in Sec. IIIII we can use a local potential 
approximation for a qualitative description of the inho- 
mogeneous density profile of the trap, by employing the 
local value of For a given value of U /i, we can there- 
fore represent the confined system by a vertical line in the 
phase diagram of the Bose Hubbard model on the square 
lattice. This representative vertical line is then shifted 
towards smaller values of t, during our constant /i scan. 
While for large values of t the system will be superfiuid, 
from Fig. Elwe expect the appearance of a central Mott 
plateau for smaller hopping amplitudes. 

In order to determine accurately the position of the 
Mott plateau formation, we monitor the evolution of the 
density in the trap center, which is shown in Fig. 1201 as a 
function of the inverse hopping, U /t. 

In accordance with Fig. ^1 the trap center has a den- 
sity larger than 1 for small values of U /t. Upon increasing 
U /t, the central density decreases, closely following the 
path of the homogeneous system. It crosses the value 1 
for U/t « 13.13, as expected from Fig. E| then under- 
goes a minimum, and reaches 1 for U/t k, 16.7, where 
a Mott plateau starts to form. Apart from the criti- 
cal region of the homogeneous 2D system, the density 



1 

0.8 
0.6 
0.4 
0.2 

-0.2 
-0.4 
-O.h 







1 

2D trap 

|Vf7=0.37 

V7t/=0.002 






superfluid 

1- - - 


trap-center 













5 








«=1 


10"; 






K//=14.3 


;7/f=io.o 


15 


n=0 








20 



t/U 



FIG. 19: Variation of fi^ in a two dimensional parabolic 
trap with curvature V/U = 0.002 for three different values of 
U/t while tuning t along a constant ^/U = 0.37 scan, shown 
within the phase diagram of the Bose Hubbard model on the 
square lattice. Along the dashed dotted line the density has 
a constant value of n = 1. 
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FIG. 20: Top: Radial density distribution of bosons within a 
two dimensional parabolic trap with curvature V/U = 0.002 
for different values of U/t for fi/U = 0.37. Bottom: Evo- 
lution of the density in the trap center while tuning t along 
a constant jj,/U = 0.37 scan. For comparison the density of 
the uniform system along the same constant /i scan is shown. 
The arrow indicates the threshold for Mott plateau forma- 
tion within the trap, which deviates by about 6% from the 
position of the Mott transition in the uniform case (vertical 
dashed line). 



in the center of the trap closely follows the behavior in 
the uniform case along the same constant fi scan. The 
rather smooth approach towards a central density of 1 is 
in agreement with the observations in Sec. IIV CI of the 
absence of critical fluctuations inside the trap. In ad- 
dition to this qualitative difference to the homogeneous 
case, the local potential approximation underestimates 
the lower value of U/t for Mott plateau formation by 
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FIG. 21: Momentum distribution function of bosons in a 
two dimensional parabolic trap with curvature V/U = 0.002 
along the line (0, Q)-{n/a, 0) in momentum space, for different 
values of U/t while tuning t along the constant fi/U = 0.37 
scan of Fig. 1201 Satellite peaks near {kx,ky) = (0.l7r/a,0) 
appear unrelated to the phase structure within the trap. 




U/t 

FIG. 22: Evolution of the coherence fraction, i.e. the height, 
n(0,0), and of the fuU width at half maximum (FWHM) 
of the coherence peak as a function of U/t while tuning t 
along the constant n/U — 0.37 scan of Fig. l20l for bosons in a 
two dimensional parabolic trap with curvature V/U = 0.002. 
The threshold for Mott plateau formation is indicated by the 
dashed line. 



about 6%, as seen from Fig.lT^ 

Having determined the value of U/ 1 for the onset of the 
central Mott plateau from measurements of the density 
distribution, we turn to a discussion of the momentum 
distribution function of bosons inside the parabolic trap 
and its evolution upon increasing the inverse hopping, 
U/t. 

In Fig. 1211 the momentum distribution function n(k) 
is shown for different values of U/t along the constant 
jj^/U = 0.37 scan of Fig. [201 From these data it is 
obvious that the presence of the fine structure peak 
near kxa/n « 0.1, reflecting the typical radial extent, 
R w 20a, of the confined Bose gas, is not related to the 
presence of a central Mott plateau, in agreement with 
earlier observations in Sec. Ivl 

Analyzing the data further, we show in Fig. |221 the 
coherence fraction, n(k = 0), and the FWHM of the co- 
herence peak as a function of U /t. In marked contrast to 
the uniform case fFig. ll7|) . but in agreement with exper- 
imental findings [l|, the coherence fraction does not dis- 
play distinct features upon emergence of a central Mott 
plateau, but instead decreases rather smoothly over a 
broader range than for the uniform system. This behav- 
ior is expected as it reflects the coexistence of both a 
Mott plateau region and a surrounding superfluid. 

Similar broadenings are observed in the evolution of 
the FWHM, which becomes rather flat in the region 
where Mott plateau formation sets in. Since changes in 
the FWHM are thus small near the threshold, care has 
to be taken when extrapolating data taken for large val- 
ues of U/t down to the flat region. For example, a linear 
extrapolation using the last two data points in Fig. |221 
would overestimate the threshold for Mott plateau for- 
mation by more than 60% of its actual value. 



However, we find that the FWHM starts to increase 
well below the threshold for Mott plateau formation. 
Furthermore, beyond this point a change in curvature 
of the graph is observed, with an inflection point located 
at the threshold point. In fact, we found the presence 
of such inflection points in the FWHM graphs at the 
transition point to be a generic feature for different trap- 
ping curvatures and dimensionalities, as will be shown 
below. This feature thus appears to be a reliable indica- 
tion for the onset of Mott plateau formation in conflned 
Bose systems. Although the FWHM is accessible exper- 
imentally 2] , limited resolution can make the location of 
the inflection point difficult, as we found the FWHM 
graphs to be rather flat in this region. 



C. Parabolic traps in 3D 

We now present results of QMC simulations of the 
Bose Hubbard model for bosons conflned in 3D parabolic 
traps. Performing an analysis like in the 2D case, we find 
that similar generic features as those obtained in 2D ap- 
ply to these 3D systems as well. 

In particular, we consider a parabolic trap with curva- 
ture V/U = 0.0125, and study the system's states for dif- 
ferent values of t/U along a line of constant fx/U = 0.25. 
For a value of U/t ~ 20, the system is still deep in the 
superfluid regime, as seen in Fig. 1231 reflecting the in- 
creased strength of the kinetic energy, due to the larger 
dimensionality. Upon decreasing the hopping t, a Mott 
plateau forms in the trap center. In Fig. [53] we trace the 
boson density in the trap center as a function of U/t, in 
order to locate the threshold for emergence of the Mott 
plateau region, indicated by the vertical arrow. Similar to 
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FIG. 23: Top: Radial density distribution of bosons within 
a three dimensional parabolic trap with curvature V/U — 
0.0125 for different values of U/t for fi/U = 0.25. Bottom: 
Evolution of the density in the trap center while tuning t 
along a constant fi/U — 0.25 scan. The arrow indicates the 
threshold for Mott plateau formation within the trap, which 
deviates by about 30% from the position of the Mott tran- 
sition in the uniform case within mean field theory (vertical 
dashed line). 



FIG. 24: Momentum distribution function of bosons in 
a three dimensional parabolic trap with curvature V/U = 
0.0125 along the line (0, 0, 0)-(7r/a, 0, 0) in momentum space, 
for different values of U/t while tuning t along the con- 
stant fi/U — 0.25 scan of Fig. 1231 Satellite peaks near 
k — (0.37r/a, 0, 0) appear unrelated to the phase structure 
within the trap. 



the 2D case, the central density approaches the value of 
1 with a flat slope. Within mean field theory '2^,'2^ and 
the local potential approximation, the threshold would 
be underestimated by about 30%, as indicated by the 
dashed vertical line in Fig. 1^ 

Analyzing the momentum distribution functions 
shown in Fig. 1241 we observe similar behavior as in the 
2D system: (i) As seen from Fig.l^S] at the threshold for 
Mott plateau formation, the coherence fraction is still 
about 10% of the overall bosonic density and decreases 
over a rather broad range of parameters. 

(ii) The FWHM of the coherence peak undergoes a 
change of curvature with an inflection point being lo- 
cated at the threshold for Mott plateau formation. The 
presence of this inflection point thus provides a robust 
indication of density restructurings also inside 3D traps, 
(iii) As in 2D, we find the presence of satellite peaks to 
be unrelated to the local density structure, as seen from 
Fig. 1^ (iv) The position of the fine structure peak is 
indicative of the spatial extent of the bosonic cloud. In 
fact, the broad peak observed at ~ 0.37r/a in Fig. 1241 
corresponds well to the radial extent i? « 6a of the Bose 
gas. 



D. Parabolic traps in ID 

Finally, we extend our analysis to the case of a ID 
parabolic trap. Fixing the chemical potential to a value of 
n/U = 0.37, similar to the 2D case, we study the system 
for different values of the hopping amplitude, t. In the 
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FIG. 25: Evolution of the coherence fraction, i.e. the height, 
n(0, 0, 0) and of the full width at half maximum {FWHM) 
of the coherence peak as a function of U/t while tuning t 
along the constant ^/U — 0.25 scan of Fig. 1231 for bosons 
in a three dimensional parabolic trap with curvature V/U = 
0.0125. The threshold for Mott plateau formation is indicated 
by the dashed line. 



upper part of Fig. 1261 the spatial density distribution is 
shown for three different values of U/t. The evolution 
of the density in the trap center as a function of U/t is 
shown in the lower part of Fig. 1261 

For U/t — 4.0 the system is in the fully superfluid 
regime, with no Mott plateaux present. Upon increas- 
ing U/t, there appears a finite regime, where two Mott 
plateaux emerge well outside the center of the trap. 
These plateaux eventually merge into an extended Mott 
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FIG. 26: Top: Spatial density distribution of bosons within a 
one dimensional parabolic trap with curvature V/U = 0.0004 
for different values of U/t, for /i/U = 0.37. Bottom: Evolu- 
tion of the density in the trap center while tuning t along a 
constant fj,/U — 0.37 scan. The threshold for Mott plateaux 
formation within the trap, and merging of the two plateaux 
are indicated by vertical arrows. 



FIG. 27: Momentum distribution function of bosons in a 
one dimensional parabolic trap with curvature V/U = 0.0004 
along the line 0-n/a in momentum space, for different values 
of U/t while tuning t along the constant ^/U — 0.37 scan of 
Fig.EHl 



plateau at a larger value of U/t. The position of these 
points is marked by the arrows in the lower part of 
Fig. 1261 This emergence of an intermediate regime with 
two well separated Mott plateaux is expected from the 
shape of the first Mott lobe in ID [2^, and the chosen 
value of /i/J7 — 0.37, and follows using a local poten- 
tial approximation, similar to the case of the 2D trap 
considered above. The reason why such an intermediate 
regime is observed in our ID simulations, but not for the 
2D case, is that in ID the largest extent of the first Mott 
lobe has a critical value of fi/U w 0.10 which is below our 
chosen value of fi/U ~ 0.37, whereas the critical value of 
IJ,/U = 0.42 in 2D is above that value. 

The corresponding momentum distribution functions 
are shown in Fig. 1271 Compared to the higher dimen- 
sional cases, the momentum distribution functions ap- 
pear broad, indicating larger incoherent contributions. 
This is expected, as in ID long range coherence cannot 
develop, even at zero temperature. Similar to the higher 
dimensional cases, we observe broad fine structure peaks 
in n{k), restricted however to U/t below the threshold 
for Mott plateau formation. At larger values of U/t, such 
fine structure is not resolved due to the large incoherent 
contribution. 

Analyzing the momentum distribution functions as 
shown in Fig. |2H1 we find that the graph of FWHM 
as a function of U/t exhibits two characteristic features: 
The increase of the slope for U/t near 5.0 corresponds 
well to the threshold for the formation of the two Mott 
plateaux. The growth of the two Mott plateaux regions 
results in the fast decrease of the coherence length in this 
regime. Beyond U/t k, 5.7, the increase in the FWHM 
is reduced, indicating that the two plateaux have merged 
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FIG. 28: Evolution of the coherence fraction, i.e. the height, 
n(0), and of the full width at half maximum {FWHM) of the 
coherence peak as a function of U /t while tuning t along the 
constant ^/U = 0.37 scan of Fig. 1261 for bosons in a one di- 
mensional parabolic trap with curvature V/U = 0.0004. The 
threshold for Mott plateaux formation within the trap, and 
for the collapse of the two plateaux into a single plateau are 
indicated by dashed lines. 



into a single plateau, which now grows at only two ends. 
The U/t dependence of the coherence fraction n(0) also 
indicates the merging of the two Mott plateaux, by a 
reduced decrease in U/t beyond U/t^ 5.7. 

Recently KoUath et al. studied the ID case us- 
ing the Density Matrix Renormalization Group method. 
Their results are in perfect agreement with our observa- 
tions. 
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VI. DISCUSSION AND CONCLUSION 

Our quantum Monte Carlo simulations provide insight 
into the physics of trapped bosonic systems on lattices. 
The validity of a local potential approximation, where lo- 
cal quantities such as the local density or compressibility 
depend mainly on the value of the local chemical poten- 
tial, is confirmed by an excellent data collapse of local 
quantities on single curves. This single curve is however 
not the same as for the homogeneous bulk system; the 
differences are particularly pronounced in the interesting 
vicinity of the transition layer between Mott insulators 
and superfluid in the parabolic trap. There the singular- 
ities due to quantum critical behavior are removed and 
replaced by smooth and broad features. 

While the behavior of the homogeneous system can 
give a qualitative overview of the phase structure realized 
locally inside the parabolic trap, quantitative results can 
only be obtained by numerical (quantum Monte Carlo) 
simulations of realistic systems. Results for realistic 2D 
parabolic traps of a size comparable to experiments have 
been presented here, as well as an analysis of the ID sit- 
uation. Three dimensional simulations have so far been 
performed only on systems with linear dimensions 2-3 
times smaller than experimental realizations, but realis- 
tic simulations will be possible in the near future using 
faster computers and improved algorithms. 

An effective ladder model, which quantitatively mod- 
els the realistic trapped system, provides clear evidence 
for the absence of quantum critical behavior in parabolic 
confinement potentials. The ladder model allows us to 
exclude the randomness imposed on the superfluid ring 
by the underlying square lattice structure as the source 
of this absence of quantum criticality. Instead, the di- 
vergences due to quantum critical fluctuations are sup- 
pressed by the inhomogeneity, and the coupling to the 
rest of the system. It will be very interesting to develop 
an effective action for this coupling, which might explain 
the power law behavior observed in Ref. JJj 21] . Fur- 
thermore, the observed absence of quantum criticality in 
both parabolic traps and ladder models with a gradient in 
the chemical potential connecting different Mott plateau 
regions calls for future investigations, and extension to 
fermionic models. 

The absence of quantum critical behavior of bosons 
in parabolic confinement potentials also agrees well with 
the fast dynamics of the "phase transition" and the ob- 
served absence of "critical slowing down" of the dynam- 
ics in experiments 0. Critical slowing down at a sec- 
ond order phase transition is caused by the long time 
scales taken to establish a uniform order parameter across 
the whole system. Small ordered domains, with differ- 
ently broken U{1) symmetry are quickly formed, but as 
these domains grow and merge, the dynamics to estab- 
lish the same U{1) symmetry breaking across neighbor- 
ing domains slows down as the domains grow in size. 
The dynamics of the "quantum phase transition" in the 
trapped system is different: driving the system from the 



superfluid to the Mott insulator happens by nucleating a 
small Mott domain inside the trap, which then grows as 
the depth of the optical lattice is increased. As the Mott 
phase grows in volume and the superfluid phase shrinks 
the "quantum phase transition" is observed. This is how- 
ever better viewed as a crossover with changing volume 
fractions of the two phases than as a phase transition: the 
large Mott plateau is always surrounded by a shell of co- 
herent superfluid. Sweeping back to the superfluid phase 
by decreasing the depth of the optical lattice, the Mott 
insulator melts and atoms join the superfluid. The dy- 
namics here is not that of two large domains merging, but 
that of a single atom joining the coherent superfluid and 
there is no critical slowing down involved in this process. 
It might be possible to experimentally observe critical 
slowing down in a trap by first driving the system deep 
into the Mott insulating region, then kicking it to destroy 
the phase coherence in the remaining superfluid shell and 
afterwards quickly driving it back into the superfluid. 

Finally, we analyzed the behavior of the momentum 
distribution function, which is accessible experimentally 
from the interference patterns of absorption images taken 
after free expansion of the atomic gas. We find that the 
full width at half maximum of the coherence peak in the 
momentum distribution function, due to its relation to 
the coherence length in the system, yields valuable infor- 
mation about density restructurings inside the trap. In 
particular, we found an inflection point in its graph upon 
increasing the lattice depth well at the threshold for cen- 
tral Mott plateau formation. Since the graph becomes 
flat in this region, detection of such features requires high 
resolution data taken in the crossover regime. 

In contrast, we found that for flat confinement poten- 
tials, realizing closed box systems, both the full width 
at half maximum, as well as the coherence fraction pro- 
vide clear signals for the Mott transition. Furthermore, 
in flat trapping geometries, quantum critical fluctuations 
are significantly more pronouned, and allow the obser- 
vation of quantum critical behavior already on optical 
lattices of currently available sizes. We thus expect the 
possible realization of flat confinement potentials to 
significantly ease the detection of true quantum critical- 
ity and the interpretation of the experimental data. 

Performing quantum Monte Carlo studies for realistic 
systems will be important for interpreting current and 
future experiments on confined Bose gases in optical lat- 
tices, and for testing our quantitative understanding of 
these systems. Such understanding will be important, 
once analog quantum computers build from fermionic 
atoms are available, for which large scale quantum Monte 
Carlo simulations will not be possible. 
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